blog
北斗三号球谐电离层用户算法代码
在BDS中全球广播电离层时延修正模型的数学结构下所示 \[{T_{ion}} = F \cdot K \cdot \left[ {{A_0} + \sum\limits_{i = 1}^9 {{\alpha _i}{A_i}} } \right]\] 电离层非播发参数形如 $\left\{ \begin{gathered} {\beta _j} = \sum\limits_{k = 0}^{12} {\left( {{a_{k,j}} \cdot \cos {\omega _k}{t_k} + {b_{k,j}} \cdot \sin {\omega _k}{t_k}} \right)} \hfill \\ {\omega _k} = \frac{{2\pi }}{{{{\text{T}}_k}}} \hfill \\ \end{gathered} \right.$ 球谐函数展开为 \[\begin{gathered} VTEC\left( {\phi ,s} \right) = \sum\limits_{l = 1}^N {\sum\limits_{m = 1}^l {{{\overline P }_{lm}}\left( {\sin \phi } \right)} } \cos \left( {ms} \right){\alpha _{lm}} + \sum\limits_{l = 1}^N {\sum\limits_{m = 1}^l {{{\overline P }_{lm}}\left( {\sin \phi } \right)} } \sin \left( {ms} \right){\beta _{lm}} \hfill \\ \quad \quad \quad \quad \quad \sum\limits_{l = 0}^N {{{\overline P }_l}\left( {\sin \phi } \right)} {\gamma _l} \hfill \\ \end{gathered} \] phi为为穿刺点的日固地磁纬度，取决于采用的坐标系。 计算球谐函数有不同的计算方法，有些方法递推误差传递较为严重，这里采用以下方法。根据 S.A. Holmes 的研究，以下方法递推过程中误差传递较弱，计算精度较高。 \[\begin{gathered} {\overline P _l}\left( u \right) = \sqrt {\frac{{2l + 1}}{{2l - 1}}} \left[ {\left( {1 - \frac{1}{l}} \right)u{{\overline P }_{l - 1}}\left( u \right)} \right. \hfill \\ \quad \quad \quad - \left. {\sqrt {\frac{{2l - 1}}{{2l - 3}}} \left( {1 - \frac{1}{l}} \right){{\overline P }_{l - 2}}\left( u \right)} \right],l \geqslant 2 \hfill \\ {\overline P _0}\left( u \right) = 1,{\overline P _1}\left( u \right) = \sqrt 3 u \hfill \\ \end{gathered} \] \[\left\{ \begin{gathered} {\overline P _{1,1}}\left( u \right) = \sqrt {3\left( {1 - {u^2}} \right)} \hfill \\ {\overline P _{l,l}}\left( u \right) = \sqrt {\frac{{2l + 1}}{{2l}}} \sqrt {1 - {u^2}} {\overline P _{l - 1,l - 1}}\left( u \right) \hfill \\ {\overline P _{l,m}}\left( u \right) = \sqrt {\frac{{\left( {2l + 1} \right)\left( {2l - 1} \right)}}{{\left( {l + m} \right)\left( {l - m} \right)}}} u{\overline P _{l - 1,m}}\left( u \right) \hfill \\ \quad \quad \quad - \sqrt {\frac{{\left( {2l + 1} \right)\left( {l - 1 + m} \right)\left( {l - 1 - m} \right)}}{{\left( {2l - 3} \right)\left( {l + m} \right)\left( {l - m} \right)}}} {\overline P _{l - 2,m}}\left( u \right) \hfill \\ l \geqslant 2,m = 1,2, \cdots ,l - 1 \hfill \\ {\overline P _{i,j}}\left( u \right) = 0,i < j \hfill \\ \end{gathered} \right.\] 由穿刺点数据计算球谐系数方法，线性方程为 \[{\mathbf{O = AX}}\] \[{\mathbf{A}} = \left[ {\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\left[ {{A_{11}}} \right]}&{\left[ {\begin{array}{*{20}{c}} {{A_{21}}}&{{A_{22}}} \end{array}} \right]}& \cdots &{\left[ {\begin{array}{*{20}{c}} {{A_{l1}}}& \cdots &{{A_{ll}}} \end{array}} \right]} \end{array}} \right]}&{\left[ {\left[ {\begin{array}{*{20}{c}} {\left[ {{B_{11}}} \right]}&{\left[ {\begin{array}{*{20}{c}} {{B_{21}}}&{{B_{22}}} \end{array}} \right]}& \cdots &{\left[ {\begin{array}{*{20}{c}} {{B_{l1}}}& \cdots &{{B_{ll}}} \end{array}} \right]} \end{array}} \right]} \right]}&{\left[ {\begin{array}{*{20}{c}} {{C_0}}&{{C_1}}& \cdots &{{C_n}} \end{array}} \right]} \\ \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots \end{array}} \right]\] \[{\mathbf{X}} = {\left[ {\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\left[ {{\alpha _{11}}} \right]}&{\left[ {\begin{array}{*{20}{c}} {{\alpha _{21}}}&{{\alpha _{22}}} \end{array}} \right]}& \cdots &{\left[ {\begin{array}{*{20}{c}} {{\alpha _{l1}}}& \cdots &{{\alpha _{ll}}} \end{array}} \right]} \end{array}} \right]}&{\left[ {\begin{array}{*{20}{c}} {\left[ {{\beta _{11}}} \right]}&{\left[ {\begin{array}{*{20}{c}} {{\beta _{21}}}&{{\beta _{22}}} \end{array}} \right]}& \cdots &{\left[ {\begin{array}{*{20}{c}} {{\beta _{l1}}}& \cdots &{{\beta _{ll}}} \end{array}} \right]} \end{array}} \right]}&{\left[ {\begin{array}{*{20}{c}} {{\gamma _0}}&{{\gamma _1}}& \cdots &{{\gamma _n}} \end{array}} \right]} \end{array}} \right]^T}\] 其中 \[\left\{ \begin{gathered} {A_{lm}} = {\overline P _{lm}}\left( {\sin \phi } \right)\cos \left( {ms} \right) \hfill \\ {B_{lm}} = {\overline P _{lm}}\left( {\sin \phi } \right)\sin \left( {ms} \right) \hfill \\ {C_l} = {\overline P _l}\left( {\sin \phi } \right) \hfill \\ \end{gathered} \right.\] 用户算法代码如下
!----------------------------------------module coment
! Version : V1.0 BDSSH 用户算法
! Coded by : song.yz
!
! Date : 2014-04-17 00:01:18 *星期五*
!-----------------------------------------------------
! Description :
! 所有程序已经得到验证
!-----------------------------------------------------
module M_inosphere
!----------------------------------------module coment
! Version : V1.0
! Coded by : song.yz
! Date : 2014-02-14 11:18:32 *星期五*
!-----------------------------------------------------
! Description :
! 地理经纬度及BDSSH表
! Post Script :
! 1.
! 2.
!
!-----------------------------------------------------
implicit real*8(a-h,o-z)
save
real*8,parameter::xlat_NP=79.85
!地磁北极地理纬度（纬度）
real*8,parameter::xlong_NP=-71.88
!地磁北极地理经度（经度）
real*8::ctable(0:24,17)
!BDSSH 周期表
end module M_inosphere
subroutine set_ctable
!---------------------------------subroutine comment
! Purpose : 设置cTable参数，目的是用于计算BDSSH中的 A0
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-02-18 10:12:08
! 设置 ctable参数
! 该函数包含在模块内部，所以直接共享模块变量
!-----------------------------------------------------
! Copyrigt (C) :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
use M_inosphere
implicit real*8(a-h,o-z)
do i=0,24
!读取系数表
read(21,*) ctable(i,:)
end do
end subroutine set_ctable
program main
!---------------------------------subroutine comment
! Purpose :
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-02-27 16:50:56
!
!-----------------------------------------------------
! Notes :
!
!
!-----------------------------------------------------
! Default File unit index :
! 11-50 input cards observation files
! 51-99 report (output) files
! 1001--9999 configuration files
!-----------------------------------------------------
! Copyrigt (C) :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
use M_inosphere
implicit real*8(a-h,o-z)
real*8::bcpara(9)
real*8,parameter::pi=3.141592653589793238462643d0
open(unit=21,file='input/a0.dat')
!a0，配置文件
!--------------------设置周期表
call set_ctable
!-----------------------------
!call test(NN,N)
xMJD=55930.04513889D0
!简约儒略日
!穿刺点地理经纬度：角度
xlat=37D00*pi/180
xlong=115.87D00*pi/180
! write(*,*) xlat,xlong,'xlat,xlong'
bcpara=(/22.685d0,-4.687d0,8.866d0, &
4.810d0, -7.183d0, -1.717d0, &
-1.337d0, 2.593d0, 1.665d0 /)
call BDSSH(cVtec,bcpara,xlat,xlong,xMJD)
write(*,*) cVtec,'vtec'
end program main
subroutine BDSSH(cVtec,bcpara,xlat,xlong,xMJD)
!---------------------------------subroutine comment
! Purpose : BDSSH模型计算理论 VTEC
! Author : Song Yezhi
!
! Versions and Changes :
! V1.0 ----2014-04-16 20:51:12
! 已经验证
!-----------------------------------------------------
! Input Parameters :
! bcpara(9)---------broadcast parameter 广播电离层参数
! xlat---------地理纬度（弧度）
! xlong--------地理经度
! xMJD---------简约儒略日
!
! Output Parameters :
! cVTEC------计算的vtec
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use M_inosphere
implicit real*8(a-h,o-z)
real*8::vrow(9),bcpara(9)
real*8::v(9) !vrow1是为了与BDSSH播发顺序匹配
NN=9
N=2
call sun_M(xlatSM,xlongSM,xMJD,xlat,xlong)
call C_A0(A0,xlatSM,xlongSM,xMJD)
call row_process_S(vrow,xlatSM,xlongSM,NN,N)
!勒让德球函数得到验证
v(1)=bcpara(1)
v(2)=bcpara(2)
v(3)=bcpara(5)
v(4)=bcpara(3)
v(5)=bcpara(4)
v(6)=bcpara(6)
v(7)=bcpara(7)
v(8)=bcpara(8)
v(9)=bcpara(9)
A1_9=dot_product(v,vrow)
cVtec=A0+A1_9
!write(*,*)vrow
!write(*,*) A0,A1_9,'A0,A1_9'
end subroutine BDSSH
subroutine row_process_S(vrow,xlatS,xlongS,NN,N)
!---------------------------------subroutine comment
! Purpose : 对一行资料进行处理，填充一行数据
! Author : song.yz
! Created : 2014-02-12 11:20:03
! Changes :
! 0. 该函数可以用来做参数估计，形成一次观测系数
! 也可以用此函数形成向量后与系数做内积，则得到
! 理论的 TVEC
!
! 1. 修改时间 2014-02-15 15:46:57
! 该程序为 row_process 的另一个版本
! row_process 处理基准是以地理经纬度为基准，并考虑
! 太阳位置的处理策略，太阳位置通过 UT得到反映
!
! row_process_M 则以日固地磁坐标系为基准，进行一行
! 数据处理，在日固地磁坐标系下，已经考虑太阳位置
! 其中太阳位置也是通过 UT 得到反映。
!
! 2. 修改时间 2014-02-15 16:12:38
! 输入单位的 xlat , xlong 单位由度 改为弧度
!
! 3. 输入的经纬度改为日固地磁坐标系，单位为弧度
! 由于坐标已经为日固地磁坐标系，所以不需要变量 UT
!
! 4. 重排vrow顺序，现在是 带谐项排完排田谐项
!
!-----------------------------------------------------
! Input Parameters :
!
! xlatM------------日固地磁纬度（弧度）
! xlongM-----------日固地磁经度（弧度）
! NN--------------参数个数
! N---------------阶数
! Output Parameters :
! vrow(NN)
!
! Subroutines used :
! *.
! *.
! *.
!-----------------------------------------------------
! Copyrigt :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
use M_inosphere
implicit real*8(a-h,o-z)
real*8,parameter::pi=3.1415926535897932384626d0
real*8::vrow(NN)
!---------------------------------------------------
real*8::PG(N,N),PZ(0:N)
u=dsin(xlatS)
!计算勒让德多项式田谐项和带谐项
call Legendre(PG,PZ,N,u)
vrow=0d0
i=1
do L=0,N
vrow(i)=PZ(L)
i=i+1
end do
do L=1,N
do M=1,L
vrow(i)=PG(L,M)*dcos(M*xlongS)
i=i+1
vrow(i)=PG(L,M)*dsin(M*xlongS)
i=i+1
end do
end do
!do L=0,N
! vrow(i)=PZ(L)
! i=i+1
!end do
end subroutine row_process_S
subroutine Rmat(A,theta,i_axis)
!---------------------------------subroutine comment
! Purpose : 坐标旋转矩阵
! Author : song.yz
! Created : 2014-02-15 13:41:46
! Changes :
! 0. 坐标系选择矩阵
! Hofmann GNSS P18
!-----------------------------------------------------
! Input Parameters :
! theta-----旋转角度 （弧度）
! iaxis-----旋转轴标号
! 1--------绕X轴旋转
! 2--------绕Y轴旋转
! 3--------绕Z轴旋转
!
! Output Parameters :
! A(3,3)
! 旋转矩阵
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8::A(3,3)
A=0d0
!设置为0
C=dcos(theta)
S=dsin(theta)
do i=1,3
A(i,i)=1d0
end do
!初始对角线设置为1
select case (i_axis)
case (1)
A(2,2) = C
A(3,2) = -S
A(2,3) = S
A(3,3) = C
case (2)
A(1,1) = C
A(3,1) = S
A(1,3) = -S
A(3,3) = C
case (3)
A(1,1) = C
A(2,1) = -S
A(1,2) = S
A(2,2) = C
case default
write(*,*) 'wrong input parameter'
end select
end subroutine Rmat
subroutine C_A0(A0,xlatS,xlongS,xMJD)
!---------------------------------subroutine comment
! Purpose :
! Author : Song Yezhi
! Versions and Changes :
! V1.0 --- 2014-02-17 14:11:27
!
!-----------------------------------------------------
! Input Parameters :
! xlatM------日固地磁坐标系纬度 （单位：弧度）
! xlongM-----日固地磁坐标系经度 （单位：弧度）
! xJD--------儒略日 （单位：天）
! Output Parameters :
! A0---------输出系数 （单位： TECu）
!
! Subroutines used :
! *. geo_legendre(PG,N,u)
! *. zone_legendre(PZ,N,u)
!-----------------------------------------------------
! Copyrigt (C) :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
use M_inosphere
!需要用到该模块的 ctable的值
implicit real*8(a-h,o-z)
real*8,parameter::pi=3.1415926535897932384626d0
!real*8::ctable(0:24,17)
!周期系数表
real*8::PER(12)
!周期，第0项的时周期为无穷大，
!在后续计算中进行考虑
real*8::beta(17),B(17)
!beta与B
real*8::PG(5,5),PZ(0:5)
!存储勒让德多项式数值，田谐项和带谐项
!周期系数，其中第0项周期为无穷大，后续计算时直接考虑
PER=(/ 1d0, 0.5d0, 0.33d0, 14.6d0, &
27d0, 121.6d0, 182.51d0, 365.25d0,&
4028.71d0, 2014.35d0, 1342.90d0, 1007.18d0 /)
!
! 系数表
! ctable
! 在模块M_inosphere中，已经在本程序启动前通过函数 set_ctable获得，并
! 保存，作为公共变量
!============================计算beta============
!xMJD=int(xMJD)+1d0/24d0
dt= xMJD-int(xMjd)
do j=0,11
tmp1=j*2d0/24d0
tmp2=(j+1d0)*2d0/24d0
if (dt>=tmp1.and.dt
! Versions and Changes :
! V1.0 ----2014-04-10 15:49:45
! 已经验证
!-----------------------------------------------------
! Input Parameters :
! xMJD-------简约儒略日
! xlat-------地理纬度 (弧度)
! xlong------地理经度 （弧度）
! Output Parameters :
! xlatSM------日固地磁坐标系下纬度
! xlongSM-----日固地磁坐标系下经度
! Subroutines used :
! *. G2M(xlat,xlong,xlatM,xlongM) 地理地磁坐标系转换
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use M_inosphere,only:xlat_NP,xlong_NP
implicit real*8(a-h,o-z)
pi=3.141592653589793238462643d0
s_lon=pi*(1-2d0*(xMJD-int(xMJD)))
xlat_N1=xlat_NP*pi/180
xlong_N1=xlong_NP*pi/180
!地磁北极经纬度化为弧度单位
tmp1=dsin(s_lon-xlong_N1)
tmp2=dsin(xlat_N1)*dcos(s_lon-xlong_N1)
call G2M(xlatM,xlongM,xlat,xlong)
xlatSM=xlatM
xlongSM=xlongM-datan2(tmp1,tmp2)
end subroutine sun_M
subroutine G2M(xlatM,xlongM,xlat,xlong)
!---------------------------------subroutine comment
! Purpose : 地理经纬度化为地磁坐标经纬度G
! Author : song.yz
! Created : 2014-02-13 17:19:36
! Notes :
! *. 本程序已经检验过
! ----------------------
!
! 地磁北极的地理经纬度在后续设置
!
! 注意本程序是把地理经纬度化为地磁经纬度
! 而没有化为日固地磁经纬度，日固地磁经纬度，在后续通过世界时
! 进行转换
!-----------------------------------------------------
! Input Parameters :
! xlat--------地理纬度（单位：弧度）
! xlong-------地理经度（单位：弧度）
!
! Output Parameters :
! xlatM-------地磁经度 （单位：弧度）
! xlongM-------地磁纬度 （单位：弧度）
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
use M_inosphere,only:xlat_NP,xlong_NP
!-----------------------------------
!xlat_NP=79.85
!地磁北极地理纬度（纬度）
!xlong_NP=-71.88
!地磁北极地理经度（经度）
!------------------------------------
implicit real*8(a-h,o-z)
real*8,parameter::pi=3.1415926535897932384626d0
real*8::XYZ(3)
!穿刺点地理坐标
real*8::XYZ_M(3)
!穿刺点地磁XYZ
real*8::Rmat1(3,3),Rmat2(3,3),Rmat3(3,3)
!旋转矩阵
xlat_N1=xlat_NP*pi/180
xlong_N1=xlong_NP*pi/180
tmpR=1d0
!单位长度，可以任意设置
XYZ(1)=tmpR*dcos(xlat)*dcos(xlong)
XYZ(2)=tmpR*dcos(xlat)*dsin(xlong)
XYZ(3)=tmpR*dsin(xlat)
!*************************************
call Rmat(Rmat1,xlong_N1,3)
call Rmat(Rmat2,pi/2d0-xlat_N1,2)
Rmat3=matmul(Rmat2,Rmat1)
XYZ_M=matmul(Rmat3,XYZ)
!**************把地理XYZ化为地磁XYZ_M
xlongM=datan2(XYZ_M(2),XYZ_M(1))
tmp=dsqrt(XYZ_M(1)*XYZ_M(1)+XYZ_M(2)*XYZ_M(2))
xlatM=datan2(XYZ_M(3),tmp)
! write(*,*) xlatM*180/pi,xlongM*180/pi,'xlatM,xlongM'
end subroutine G2M
subroutine Legendre(PG,PZ,N,u)
!---------------------------------subroutine comment
! Purpose : 递推法计算勒让德多项式带谐项和田谐项
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-02-17 15:11:24
! 递推法计算正交归一化勒让德多项式
! 算法详见《航天器轨道理论》
! 田谐项和带谐项解耦
! V2.0 ----2014-02-27 17:10:28
! 把田谐项和带谐项写为内部函数
! 语法参考《Fortran 95/2003 for scientists and Engineers》 Chapman
! chapter 9.7 Internal Procedures Page 446
! 递推过程参考V1.0公式，与以下公式同
! S.A. Holmes Jounal of Geodesy 2002
! A unitfied approach to the Clenshaw summation and recursive
! computation of very high degree and order normalised associated
! Legendre functions
!-----------------------------------------------------
! Input Parameters :
! N-----------阶数
! u-----------自变量
! 注意：
! u----------在势函数理论中（如牛顿引力势、静电场势等）
! u-实际为 sin(phi)
! phi为展开点在球坐标下的纬度
! Output Parameters :
! PG(N,N) -------归一化勒让德多项式田谐项
! PZ(0:N)---------归一化勒让德多项式带谐项
! 注意：
! PZ的指标从0阶开始，0阶为常数项
! PG的指标从1阶开始
! Subroutines used :
! *. geo_legendre(PG,N,u)
! *. zone_legendre(PZ,N,u)
!-----------------------------------------------------
! Copyrigt (C) :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8::PG(N,N),PZ(0:N)
call geo_legendre(PG,N,u)
call zone_legendre(PZ,N,u)
contains
subroutine geo_legendre(PG,N,u)
!---------------------------------subroutine comment
! Purpose : 计算归一化勒让德多项式田谐项
! Author : song.yz
! Created : 2014-02-10 20:06:19
! Changes :
! *.
!-----------------------------------------------------
! Input Parameters :
! N------阶数
! u------自变量
! Output Parameters :
! PG（N,N） 归一化勒让德多项式
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8::PG(N,N)
PG=0D0
PG(1,1)=dsqrt(3*(1-u*u))
!计算所有扇谐项
do L=2,N
TMP=(2d0*L+1)/2D0/L*(1-u*u)
TMP=dsqrt(TMP)
PG(L,L)=TMP*PG(L-1,L-1)
end do
!从第二列开始计算，至倒数第二列，最后一列为扇谐项已经计算完毕
!为了提高计算效率，第一列暂不计算，对角线以上的元素为0，前面已经设置完毕
!第一列要用数组以外的0元素
do M=2,N-1
do L=M+1,N
tmp1=(2d0*L+1D0)*(2D0*L-1D0)/(L+M)/(L-M)
tmp1=dsqrt(tmp1)*u
tmp2=(2d0*L+1)*(L-1D0+M)*(L-1D0-M)/(2D0*L-3D0)/(L+M)/(L-M)
tmp2=dsqrt(tmp2)
PG(L,M)=tmp1*PG(L-1,M)-tmp2*PG(L-2,M)
end do
end do
! 计算第一列
M=1
do L=2,N
tmp1=(2d0*L+1D0)*(2D0*L-1D0)/(L+M)/(L-M)
tmp1=dsqrt(tmp1)*u
tmp2=(2d0*L+1)*(L-1D0+M)*(L-1D0-M)/(2D0*L-3D0)/(L+M)/(L-M)
tmp2=dsqrt(tmp2)
if (L==2) then
PG(L,M)=tmp1*PG(L-1,M)
else
PG(L,M)=tmp1*PG(L-1,M)-tmp2*PG(L-2,M)
end if
end do
end subroutine geo_legendre
subroutine zone_legendre(PZ,N,u)
!---------------------------------subroutine comment
! Copyright : Shanghai Astronomical Observatory
! Author : song.yz
! Purpose :
! Version & Changes :
! V 1.0- (Created)---------- 2014-01-14 15:31:46
! 归一化Legendre带谐项多项式计算函数
! 《航天器轨道理论》---P109
! 注意：带谐项从0开始，这样可以直接与球谐展开中的常数项合并
! 而球谐项则从1开始，如果球谐项从0的话，则其中有一项是冗余的
!-----------------------------------------------------
! Input Parameters :
! n-------------最大带谐阶数
! u-------------自变量
!
! Output Parameters :
! PZ(0:N)----------归一化legendre带谐项函数值
!
! Subroutines used :
! *.
! *.
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8::PZ(0:N)
PZ(0)=1D0
PZ(1)=dsqrt(3d0)*u
do i=2,N
TMP1=(2D0-1D0/i)*u
TMP2=dsqrt((2D0*i-1D0)/(2D0*i-3D0))*(1D0-1d0/i)
TMP3=dsqrt((2D0*i+1)/(2d0*i-1))
PZ(i)=TMP3*(TMP1*PZ(i-1)-TMP2*PZ(i-2))
end do
end subroutine zone_legendre
end subroutine Legendre
subroutine cal2MJD (I,M,K,H,TJD,TMJD)
!---------------------------------subroutine comment
! Purpose : computes JD and MJD from calendar date and time
! Author : Novas .
! Modified by : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-12 09:43:03
!-----------------------------------------------------
! Input Parameters :
! I = YEAR (IN)
! M = MONTH NUMBER (IN)
! K = DAY OF MONTH (IN)
! H = UT HOURS (IN)
!
! Output Parameters :
! TJD = JULIAN DATE (OUT)
! TMJD = MJD (out)
!
! Ref :
! *. Novas ----- JULDAT
! *.
! 如
! Year M D MJD
! 2008 01 01 54466
! 2008 01 02 54467
! 2008 01 03 54468
!----------------------------------------------------
real*8:: H,TJD,TMJD
! JD=JULIAN DAY NO FOR DAY BEGINNING AT GREENWICH NOON ON GIVEN DATE
JD = K-32075+1461*(I+4800+(M-14)/12)/4+367*(M-2-(M-14)/12*12)/12 &
-3*((I+4900+(M-14)/12)/100)/4
TJD = JD - 0.5D0 + H/24.D0
TMJD=TJD-2400000.5D0
end subroutine cal2MJD